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1  Introduction 


As  part  of  The  National  Environmental  Policy  Act  of  1969,  it  is  “the  continuing 
responsibility  of  the  Federal  Government  to  use  all  practicable  means”  to  study,  analyze  and 
report  on  the  impact  military  operations  have  on  the  environment.  This  Environmental 
Impact  Analysis  Process  (EIAP)  is  managed  at  the  Air  Force  Space  and  Missile  Systems 
Center  (SMC)  by  the  Acquisition  Civil  and  Environmental  Engineering  Group  (SMC/AXF.) 
For  more  than  the  past  ten  years,  SMC/AXF  has  conducted  a  multi-phase,  multi-year  study  of 
the  effects  of  space  launch  vehicle  (and  object  reentry)  sonic  booms  have  on  land  and  sea 
mammals.  This  requires  both  computer  models  that  predict  the  in-air  sonic  boom  that 
reaches  the  Earth’s  surface  from  launch  and  reentry  vehicles,  as  well  as  computer  models  that 
predict  sonic  boom  pressure  disturbances  that  are  expected  under  the  ocean’s  surface. 

As  part  of  the  SMC/AXF  sonic  boom  study,  a  “sonic  boom  wavy  surface”  computer  code 
was  developed  by  HKC  Research  at  USC  over  the  last  few  years  and  was  distributed  in  April 
2003  via  CD  (Ref.  1)  through  PARSONS,  Pasadena,  to  SMC/AXF.  A  user’s  guide  (Ref.  2) 
for  the  code  was  provided  to  demonstrate  the  steps  of  execution  required  to  obtain  a  solution 
for  a  given  set  of  input  parameters.  The  theory  upon  which  this  code  is  based  is  presented  in 
Ref.  3.  Examples  of  and  extensions  to  the  theory  are  given  in  Ref.  4.  A  final  report  was  also 
distributed  through  PARSONS,  Pasadena  (Ref.  5.) 

The  purpose  of  the  current  report  is  to  document  the  code,  using  notation  which  matches  (as 
closely  as  possible)  that  which  is  used  in  the  theoretical  development  put  forth  in  Ref.  3. 

This  work  thus  provides  a  correspondence  between  Cheng  and  Lee’s  equations  and  the 
routines  which  are  executed  within  the  code.  In  addition,  code  results  are  plotted  for  a 
sample  input  set  in  order  to  provide  benchmark  results  for  future  users  and  for  possible  future 
releases  by  Cheng  and  Lee. 

The  emphasis  of  this  report  is  on  code  documentation.  No  assessment  nor  interpretation  of 
the  theory  is  offered.  Readers  interested  in  the  theory  may  consult  Ref.  3. 

This  report  is  intended  to  supplement  the  user’s  guide  (Ref.  2)  supplied  by  Cheng  and  Lee  in 
that  it  provides  a  bridge  between  the  theoretical  equations  in  Ref.  3  and  the  code  itself  (Ref. 
1.)  The  user’s  guide  (Ref.  2)  supplied  on  the  CD  does  not  provide  this  bridge.  HKC 
Research  has  reviewed  our  work  and  approved  of  this  supplemental  documentation. 

The  Cheng  and  Lee  code  is  not  currently  used  as  a  standard  tool  for  predicting  sonic  boom 
transmission  into  the  ocean.  The  more  standard  approach  is  due  to  Sawyers  (Ref.  6),  wherein 
the  effect  of  a  wavy  surface  is  neglected. 
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2  Problem  Description 


The  problem  attacked  in  Ref.  3  is  that  of  a  sonic  boom,  generated  in  air,  impinging  upon  and 
interacting  with  a  water  surface.  It  is  desired  to  find  the  overpressure  (as  a  function  of  time) 
at  some  specified  (fixed)  location  beneath  the  water’s  surface.  For  the  case  of  a  flat  water 
surface  (no  waves),  this  problem  was  worked  out  some  time  ago  by  Sawyers  (see  Ref.  6.) 

The  generalization  considered  in  Ref.  3  is  the  situation  where  there  is  a  wave  train  on  the 
water  surface.  The  water  surface  waves  are  assumed  to  propagate  in  a  direction  which  is 
nearly  opposite  to  the  horizontal  flight  path  of  the  supersonic  body  which  generates  the  sonic 
boom.  If  the  supersonic  body  speed  is  denoted  by  U ,  the  water  wave  (phase)  speed  by  c,  and 
the  distance  between  water  wave  crests  by  A,  ,  then  it  is  easy  to  see  that  the  sonic  boom  “hits” 
the  water  wave  crests  with  a  frequency  (in  hz)  of  (f/  +  c)/2.  And  since  the  water  wave 
phase  speed  is  very  small  compared  with  U,  the  interaction  frequency  is  closely  approximated 
hy  U  !  A  .  This  interaction  frequency  is  an  important  parameter  in  the  analysis  contained  in 
Ref.  3. 

In  the  Cheng  &  Lee  theory,  the  desired  overpressure  time  history  is  assumed  to  be  comprised 
of  two  contributions. 


pit,x,z)  =  p^{t,x,z)  +  p2(t,x,z) 


(2.1) 


[In  the  above  equation,  over-bars  are  used  to  denote  dimensional  variables. 

The  reason  for  this  notation  is  that  in  the  following  analysis,  we  will  almost 
always  deal  with  dimensionless  (scaled)  variables,  which  will  be  denoted 
without  the  over-bar.] 

The  symbol  p  is  used  to  represent  the  overpressure  (psf),  as  a  function  of  time  T  (sec),  at  a 
fixed  depth  of  7  (ft)  and  a  fixed  horizontal  location,  x  (ft).  The  first  term,  p^ ,  is  the  flat- 
surface  (Sawyers)  contribution,  and  the  second,  ,  is  the  contribution  due  to  interaction  of 
the  sonic  boom  with  the  wavy  water  surface.  As  shown  in  Ref.  3,  the  flat  surface  term 
dominates  at  small  depths.  But  at  larger  depths,  the  interaction  term,  p^ ,  dominates.  This  is 
due  to  the  different  dependence  of  the  two  terms  with  depth,  when  depth  is  large. 
Speeifically,  these  two  terms  fall  off  with  depth  as  follows: 

Pi  z  whereas  P2^z 


The  geometry  of  the  problem  is  explained  with  the  aid  of  Fig.  I,  which  shows  a  snapshot  at 
F  =  0 .  At  this  instant,  the  surface  signature  of  the  sonic  boom  (which  is  running  to  -x  at 
speed  1/ )  is  non-zero  only  on  the  interval  0<x<L. 
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2.1  Nomenclature 


The  following  nomenclature  is  used  for  dimensional  variables  of  interest. 

F  =  time  (sec) 

Z  =  depth  (ft) 

^  =  horizontal  position,  measured  positive  opposite  to  flight  direction  (ft) 

p  =  overpressure  (psf) 

^max  =  maximum  overpressure  in  sonic  boom  signature  at  the  water  surface  (psf) 

U  =  horizontal  flight  speed  of  supersonic  body  (ft/sec) 

L  =  sonic  boom  signature  length  at  water  surface  (ft) 

=  speed  of  sound  in  air  (ft/sec)  [Code  assumes  =1118.2  ft/sec.] 
r  =  ratio  of  water  sound  speed  to  air  sound  speed  [Code  assumes  r  =  4.53] 

A  =  distance  between  water  wave  crests  (ft) 

k  =  Itt  I A  =  water  wave  number  (ft'*) 

\  =  water  wave  crest  height  (ft) 

c  =  water  wave  phase  speed  (ft/sec) 

A  =  swept  angle  of  impact  point  (off-track  location)  (deg) 

'F  =  non-alignment  angle  between  water  wave  propagation  and  flight  direction  (deg) 

Q,  =  interaction  frequency  of  sonic  boom  with  crests  (rad/sec) 

[The  lengths,  L,  A,  and  \  ,  and  the  angles,  A  and  T ,  are  explained  via  Figs.  1  and  2.] 

Some  corresponding  non-dimensional  variables  are  defined  as  follows. 

t  =  T  U I L  (dimensionless  time) 

Z  -TIL  (dimensionless  depth) 

X  =  X  I L  (dimensionless  horizontal  position) 

p  =  P  I  -fJnax  (dimensionless  overpressure) 

=  U I  (Mach  number  on  air  side  of  surface) 
k  =  k  L  (dimensionless  wave  number) 

S  =  A^l  A  (slope  parameter) 

Q.  =  Q.  LIU  (dimensionless  interaction  frequency) 


4 


2.2  Model  Parameters 


Various  parameters  which  appear  in  the  Cheng  &  Lee  model  are  defined  in  the  table  below. 
The  first  column  contains  an  equation  number  pertaining  to  the  current  document.  The  third 
and  fourth  columns  contain  the  cross  reference  to  Ref.  3. 


Eq.  No. 

Definition 

Page 
(Ref.  3) 

Equation 

(Ref.  3) 

D.l 

=  IJ la (assume  =1 1 18.2  ft/sec) 

5 

line  11 

D.2 

=M^cosA 

13 

above  (5.3) 

D.3 

=  ^cos(A  +  'F) 

13 

(5.4b) 

D.4 

it2=;fcsin(A  +  'F) 

13 

(5.4b) 

D.5 

13 

below  (5.3) 

D.6 

P  =  2i!:M„^cos(A  +  'T) 

14 

(5.7e) 

D.7 

14 

(5.7f) 

D.8 

a  =  4P^-ABlQl[2Bl) 

14 

(5.8b) 

D.9 

M  =  Pl(2Bla) 

15 

(5.10d) 

D.IO 

Mjy  =M^lr  ( assume  r  =  4.53  ) 

5 

line  10 

D.ll 

11 

(4.1c) 

D.12 

P^=lkMl  cos(A  +  'T) 

D.13 

(2w  =^'[M^-(l  +  M^)sin'(A  +  'F)] 

D.14 

+45^  Qw 

19 

(5.20b) 

D.15 

[Equivalent  to  Cheng:  z-^„z] 

21 

above  (5.26) 

5 


2.3  Interaction  Frequency 

The  interaction  frequency  is  given  in  terms  of  wave  number  by  the  following  relation  (see 
Ref.  3,  Eq.  5.2b.) 


Q.  =  ck  ■>rU 

Or,  using  definition  D.3  in  the  above  table  (see  Ref.  3,  Eq  5.4b), 

Q  =  it[c  +  f/cos(A  +  'P)] 


(2.2) 


(2.3) 


In  non-dimensional  form,  this  relation  is 

Q  =  k[c/U+cos(A-\-^)]  (2.4) 

Within  the  code  (and  most  places  within  the  theoretical  development  in  Ref.  3)  the  first  term 
in  the  brackets  is  neglected,  since  the  water  surface  wave  phase  speed  is  much  smaller  than 
the  boom  speed  in  air.  Thus  the  frequency  relation  actually  used  in  the  code  is 

Q  =  itcos(A  +  'P)  =  )fe,  (2.5) 
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3  Cheng  &  Lee  Theory 


With  the  preliminaries  of  the  first  two  sections,  the  theory  of  Ref.  3  will  now  be  summarized 
in  a  manner  which  places  it  in  correspondence  to  the  code  contained  in  Ref.  1. 

The  code  computes  the  two  overpressure  contributions  in  the  following  non-dimensional 
form. 


p{x,z)  =  Pi{x,z)  +  P2{x,z) 


(3.1) 


3.1  Flat-surface  Term 

The  first  term  is  the  flat-surface  contribution,  which  is  calculated  according  to  the  Sawyers 
theory  (see  Ref.  6)  as 


Py{x,z)=  -u(x,z) 


(3.2) 


where  the  function  u(x,  z)  is  computed  within  the  code  as 


u{x,z)  =  Re 


-  j  -Mo(j:)[Log(Zo  -1)-Log(zo)]| 

zr  ZQ-g  I 


(3.3) 


where  “Re”  means  “real  part”  and  Zq  is  the  complex  position: 


Zo  =  X  -I-  i  Zc  (where  Zs=B„z) 


(3.4) 


The  function,  Uq{x)  ,  is  the  unit  waveform  of  the  sonic  boom  at  the  water  surface.  This 

function  is  non-zero  only  for  0  <  ;c  <  1 .  Two  “canned”  surface  waveforms  are  provided 
within  the  code:  one  for  an  N-wave  (see  Fig.  3)  and  another  for  a  Titan  IV  (see  Fig.  4.)  In 
addition,  the  user  may  alter  the  routine  for  Uq(x)  ,  and  thereby  specify  any  arbitrary  unit 
waveform  as  input. 
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3.2  Wavy  Surface  Term 

The  wavy  surface  term  is  expressed  as  follows.  [See  Ref.  3,  Eq.  (5.11).] 

P2(x,z)  =  ^  (3-5) 

In  the  above  equation,  S  is  the  slope  parameter  defined  in  Section  2.1.  In  the  code,  the  phase 
factor  is  computed  using  non-dimensional  variables  as  follows. 

P2(x,z)^SRe^  P2  (3.6) 


3.2.1  Complex  Potential  at  Surface,  ^2^ 

Calculation  of  the  complex  overpressure  amplitude,  P2  ,  at  a  field  location  ( x,  z )  involves 
several  steps.  One  begins  by  computing  the  surface  ( z  =  0)  values  for  the  complex  potential 
02  associated  with  this  overpressure.  Specifically,  the  code  computes 

02(x)  =  2^Six)  (3.7) 


where 


Six)  =  ^H(x^)G{x-x^)dx^  for  0<x<l  (3.8) 

0 

and 

1 

S(a:)  =  j//(x,)G(x-x,)d[x:,  for  x>l  (3.9) 

0 

where  the  integrand  is  defined  by  the  product  of  two  functions: 

H{x{)^  ik^u^ (x, )  -  BIuq (x,  )  (3.10) 

where  is  the  derivative  of  the  unit  surface  waveform,  and 

G(^)  =  e-'<*--^“>Oo(«^)  (3.11) 

The  latter  function  (see  Ref.  3,  Eq.  5.10c)  involves  7q  ,  the  zero-order  Bessel  function  of  the 
first  kind.  It  can  be  shown  that  Eq.  (3.7)  [along  with  Eqs.  (3.8,  10,  and  11)]  is  exactly 
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equivalent  to  Eq.  (5.8a)  on  page  14  of  Ref.  3,  if  in  the  latter  equation,  one  sets  z  =  0  and 
identifies  f'  as  -u^. 

3.2.2  Interaction  Overpressure  at  Surface, 

At  the  surface,  the  overpressure  due  to  interaction  is  computed  by  a  surface  routine  in  the 
code  as  follows. 


KsurfW  =  - 


dx 


The  above  expression  corresponds  to  Eqs.  (5.8e,f)  in  Ref.  3. 


(3.12) 


NOTE:  The  surface  routine  provided  on  the  CD  contains  two  minor  errors.  It  reverses  the 
sign  in  (3.12)  and  incorrectly  uses  k  rather  than  in  the  computation.  The  output 
for  the  sample  problem  given  in  Section  5  incorporates  the  sign  correction. 


3.2.3  Complex  Transformed  Amplitude, 

A 

Having  computed  ^2(^)  >  ^  complex  function  is  now  computed  which  is  part  of  the  Fourier 
transform  of  the  desired  function,  P2ix,  z) .  Within  the  code,  this  function  is  computed  as 
follows. 


A(a  =  A/^)+A,(^) 


(3.13) 


where 


and 


\  y 


0 


As(i)  =  2 


y 


4^ 


(3.14) 


(3.15) 


where,  in  (3.15),  the  derivative  of  the  Sawyers  solution  with  depth,  evaluated  at  the  surface, 
has  been  denoted  by 


du. 

qix)  s  — I 
dz 


z=0 


(3.16) 
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The  computation  of  A(^)  via  Eqs.  (3.13)  -  (3.16)  is  equivalent  to  Eqs.  (5.18a,b)  on  page  18 
of  Ref.  3. 

3.2.4  Underwater  Interaction  Overpressure. 

With  the  transformed  amplitude  A(^)  obtained  in  section  3.2.3,  the  overpressure  due  to 
wavy  surface  interaction  is  computed  as  the  following  inverse  Fourier  transform.  [See  Ref.  3, 
Eq.  (5.17).] 


P2  z)  =  -4=  I A(^)  cr{^,  z)d^ 
s27r  ji 


(3.17) 


where 


exp 


exp 


z] 


for  Ki<^)<0 
for  K(^)>0 


and 


(3.18) 


(3.19) 


and  ,  and  Qiv  ^re  defined  in  Eqs.  D.  1 1  -  D.13  in  Section  2.2.  Note  that  Eqs.  (3.18) 

and  (3.19)  above  correspond,  respectively,  to  Eqs.  (5.15b)  and  (5.13b)  on  page  17  of  Ref.  3. 


3.2.5  Far-Field  (Deep  Water)  Solution.  ^2,farU-z) 

In  addition  to  the  exact  solution  computed  as  the  inverse  Fourier  transform  in  Eq.  (3.17),  the 
code  also  computes  an  asymptotic  approximation,  which  is  valid  for  large  z.  This  is 
computed  as 


F2.farUz) 


V^A(^>) 


exp(/£’) 


(3.20) 


where 


E  = 


-’w 


y2B^  J 


7t 

Zc - - 


(3.21) 


1] 


Eqs.  (3.20)  and  (3.21)  are  equivalent  to  Eq.  (5.26)  on  page  21  of  Ref.  3.  The  parameters 
and  Zs  in  the  above  equations  are  defined  in  Eqs  D.14  and  D.15  in  Section  2.2.  In  Eqs. 

(3.20)  and  (3.21),  two  additional  definitions  are  required.  These  are  given  below,  and 
correspond  to  Eqs.  (5.23c)  and  (5.24)  on  page  21  of  Ref.  3. 

t]  =  xl(By^z)  =  xl  Zs  (3.22) 


P^-S^T}/yll  +  TJ^ 

2Bl 


(3.23) 
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4  Cheng  &  Lee  Code 


4.1  Contents  of  the  Distribution  CD 

The  code  which  performs  the  calculations  described  in  the  previous  section  was  received  on  a 
CD  (Ref.  1.)  The  relevant  files  on  that  CD  are  discussed  below. 

4.1.1  FORTRAN  Source  Code 

There  are  six  files  on  the  CD  which  contain  the  FORTRAN  source  code  required  to  perform 
the  calculations.  (See  the  following  table.) 


Step 

File 

Purpose 

See 

Section: 

1 

2d_gen26.for 

Computes  the  Sawyers  solution,  u{x,  z) 

3.1 

2 

Sw_stl6.for 

Computes  surface  potential,  ^2^ 

3.2.1 

3 

SUR_P26.FOR 

Computes  interaction  overpressure  at  surface,  P2,surf 

3.2.2 

4 

Sw_st27.for 

Computes  amplitude,  A(^) 

3.2.3 

5 

Sw_st34.for 

Computes  P2{x,z)  and  P2,iAx,z) 

3.2.4  & 

3.2.5 

df2.for 

Function  which  defines  unit  surface  waveform,  Uq{x) 

3.1 

The  above  six  files  were  copied  to  a  folder,  and  the  latter  was  renamed  to  “df.for”  in  order  to 
be  recognized  in  the  “include”  statement  in  the  other  modules.  The  first  five  programs  in  the 
above  table  were  then  compiled  (under  Lahey  95.)  These  five  codes  were  then  run 
sequentially,  using  inputs  to  the  screen  for  the  sample  problem  described  in  the  User’s  Guide 
(Ref.  2.)  The  text  files  output  by  these  five  codes  were  examined.  [NOTE:  For  the  sample 
problem,  the  selected  unit  surface  waveform,  Uq{x)  ,  is  an  N-wave.] 

4.1.2  Supplied  Executables 

The  CD  also  contains  ten  executable  files.  However,  only  five  of  them  are  relevant.  The 
others  apparently  represent  an  earlier  version  of  the  code.  The  relevant  ones  are  given  in  the 
second  column  of  the  following  table. 


Step 

i.-WBWIffilWI.Iil 

Corresponds  to  Source  File: 

Test  Results  Agree? 

1 

2d_gen26.for 

Yes 

2 

Sw_stl5.exe 

Sw_stl6.for 

Yes 

3 

SUR_P26.exe 

SUR_P26.FOR 

Yes 

4 

Sw_st27.for 

Yes 

5 

Sw_st34.for 

Yes 

In  the  above  table,  the  last  column  means  that  the  executable  supplied  on  the  CD  was  run 
with  the  same  inputs  as  were  used  for  the  executable  generated  by  compiling  the  source  file. 
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A  “Yes”  entry  means  the  two  executables  produced  the  same  results.  It  is  clear  from  these 
tests,  that  the  six  FORTRAN  files  contain  the  source  code  which  was  used  to  generate  the 
five  relevant  executables  (second  column  of  above  table)  supplied  on  the  CD. 

4.1.3  EXCEL  Spreadsheets 

There  are  three  relevant  EXCEL  files  on  the  CD.  These  are  used  to  process  and  plot  output 
from  the  five  execution  steps  given  above.  The  following  table  gives  the  utility  of  the 
spreadsheets. 


EXCEL  Spreadsheet 

Requires  Pasting  Output  From: 

Plot_Phi_Fta.xls 

phi_st.txt  (Step  2)  and  fta_st.txt  (Step  4) 

p2_zeq0.txt  (Step  3)  and  p0_zeq0.txt  (Step  3) 

P_Zeq_l.xls 

p0_zl.txt  (Step  1) ,  p2_st.txt  (Step  5),  and  p2_far.txt  (Step  5) 

4.1.4  User’s  Guide 


The  final  file  of  interest  on  the  distribution  CD  is  “User_Guide_IV.ppt”  which  contains  the 
User’s  Guide  (Ref.  2.) 

4.2  Running  a  Sample  Problem 

The  steps  required  to  obtain  the  solution  to  a  sample  problem  are  now  described.  Screen 
inputs  are  given  to  the  right  of  the  equal  signs  (courier  font).  [All  numerical  results  which 
follow  were  generated  using  compiled  source  code.  However,  the  resulting  executables  have 
been  shown  in  Section  4.1.2  to  be  equivalent  to  the  executables  supplied  on  the  CD.] 

4.2.1  Step  1;  Generate  the  Sawyers  Solution 

Run  the  code  for  the  Sawyers  solution,  “2d_gen26.exe.” 

Screen  Inputs: 

output  file  name,  = 

Mach  number  at  the  air  side,  = 

Begin  calculation  at  Xbegin  = 

End  calculation  at  Xend  = 

Number  of  divisions,  Nx  = 

Scaled  Depth,  z  = 


pO_zl . txt 
1.5 
-20.0 
20.0 
200 
0.5 


Output: 

CD  Two-column  text  file  “p0_zl.txf’  which  contains: 

X  M(r,z  =  0.5) 

from  X  =  Xbegin  to  jc  =  Xend  at  intervals  of  Ar  =  (Xend  -  Xbegin)/Nx  . 
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4.2.2  Step  2;  Compute  Potential  at  Surface 

Run  the  code  for  the  surface  potential,  “Sw_stl6.exe.” 

Screen  Inputs: 

Wave  number,  k  = 

Mach  number  at  the  air  side,  = 

Swept  angle  (deg) ,  A  = 

Non-alignment  angle  (deg) ,  'F  = 

End  of  Phi  calculation,  Xstop  = 
Integration  step,  Xstep  = 


14.0 

1.5 

0.0 

0.0 

20.0 

0.01 


Output: 


Q  Three-column  text  file  “phi_st.txt”  which  contains  real  and  imaginary  parts  of 


f“  A  “I  f~  ^ 

Re  02^^)  I™  02^^) 


from  r  =  0  to  jc  =  Xstop  at  intervals  of  A  jc  =  Xstep . 


CD  Unformatted  (binary)  file  “phi_dat”  needed  for  steps  3  and  4. 

4.2.3  Step  3;  Compute  Interaction  Overpressure  at  Surface 

Run  the  code  for  the  surface  overpressure,  “SUR_P26.exe.” 

Screen  Inputs: 

input  file  name,  =  phi_dat 


Output: 

CD  Two-column  text  file  “p0_zeq0.txt”  which  contains: 

X  u(x,z  =  0)  =  Uq(x) 

from  r  =  0  to  X  =  Xstop  at  intervals  of  A  x  =  5*Xstep . 

Q  Three-column  text  file  “p2_zeq0.txt”  which  contains  real  and  imaginary  parts  of 
KsurfW: 

^  Re[p2^^(x)]  Im[p2_,„rf(x)] 

from  X  =  0  to  X  =  Xstop  at  intervals  of  A  x  =  5*Xstep . 
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4.2.4  Step  4;  Compute  Transformed  Amplitude 

Run  the  code  for  the  amplitude,  “Sw_st27.exe.” 

Screen  Inputs: 

input  Phi  file  name. 

Begin  calculation  at  Xibegin 
End  calculation  at  Xiend 
Integration  step,  dxi 
End  of  Phi  Calculation,  Xstop 


=  phi_dat 
=  -20.0 
=  20.0 
=  0.01 
=  20.0 


Output: 

CH  Three-column  text  file  “fta_st.txt”  which  contains  real  and  imaginary  parts  of  A(^) : 

^  Re[A(^)]  Im[A(^)] 

from  ^  =  Xibegin  to  ^  =  Xiend  at  intervals  of  A  ^  =  dxi . 

Cll  Unformatted  (binary)  file  “fta_dat”  needed  for  step  5. 

4.2.5  Step  5;  Compute  Interaction  Overpressure  at  Depth 

Run  the  code  to  compute  the  wavy  surface  overpressure  at  depth,  “Sw_st34.exe.” 

Screen  Inputs: 

Input  A(xi)  file  name. 

Scaled  Depth,  z 
Begin  calculation  at  Xbegin 
End  calculation  at  Xend 
Number  of  divisions,  Nx 

Output: 

CD  Three-column  text  file  “p2_st.txt”  which  contains  real  and  imaginary  parts  of  P2{x, z) 
X  Re[p2(^,t)]  Im[p2(jr,z)] 
from  r  =  Xbegin  to  r  =  Xend  at  intervals  of  Ar  =  (Xend  -  Xbegin)/Nx  . 

CD  Three-column  text  file  “p2_far.txt”  which  contains  real  and  imaginary  parts  of 
P2m{x,z)-. 

^  Re[p2f^(r,z)]  Im[p2f„(r,z)] 

from  x  =  Xbegin  to  x  =  Xend  at  intervals  of  A  r  =  (Xend  -  Xbegin)/Nx  . 


=  fta_dat 
=  0.5 

=  -20.0 
=  20.0 
=  200 
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4.2.6  Step  6;  Examine  Output  Usine  Spreadsheets 

Examine  the  output  of  Steps  1-5  using  the  spreadsheets  provided. 

Spreadsheet  1:  “Plot_Phi_Fta.xls” 

Paste  the  output  in  the  files  “phi_st.txt”  (Step  2)  and  “fta_st.txt”  (Step  4)  into  this 
spreadsheet  to  look  at  plots  of  the  real  and  imaginary  parts  of  ^2(^)  A{g) .  This 
spreadsheet  performs  no  calculations. 

Spreadsheet  2:  “P_Zeq_l.xls” 

This  is  the  most  important  spreadsheet.  It  allows  display  of  the  final  result,  p(f,z) ,  given  in 
Eq.  (2.1).  (The  horizontal  field  location  is  always  taken  as  J  =  0 .)  To  use  this  spreadsheet, 
paste  the  three-column  output  in  “p2_st.txt”  (Step  5)  into  columns  B,  C,  and  D.  Paste  the 
three-column  output  in  “p2_far.txt”  (Step  5)  into  columns  F,  G,  and  H.  Finally,  paste  the 
two-column  output  in  “pO_zl.txt”  (Step  1)  into  columns  J  and  K.  At  the  top  of  the 
spreadsheet,  it  is  also  necessary  to  enter  the  values  of  the  parameters: 
and  L. 

In  cell  K3,  the  non-dimensional  interaction  frequency  is  computed  according  to  Eq.  (2.5). 

Q.  =  k  cos(A  +  'P)  = 

In  cell  K5,  the  velocity  of  the  supersonic  body  is  computed: 

U  =  a^M^  (assumption:  1118.2  ft/sec) 

In  column  M,  the  spreadsheet  transforms  the  non-dimensional  horizontal  position,  x,  in 
column  B  to  dimensional  time  as  follows. 

T  ^xL/U 


In  column  N,  the  flat  surface  overpressure  is  computed  from  the  Sawyers  result  in  column  K. 
[SeeEq.  (3.2).] 


p,{t  ,z)=-P^M^,z) 

In  column  O,  the  wavy  surface  overpressure  is  computed  from  the  position,  jc,  in  column  B 
and  the  real  and  imaginary  parts  of  P2{x, z)  in  columns  C  and  D.  [See  Eq.  (3.6).] 

P2(t,z)  =  S  Re[p2(:f,z)e''’“''] 
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In  column  P,  the  final  result,  i.e.,  the  total  overpressure  at  depth,  is  tabulated  as  the  sum  of 
column  N  and  O.  [See  Eq.  (2.1).] 

p(t^z)  =  Px(t,z)  +  P2(J,z) 

Spreadsheet  3:  “P_Zeq_0.xls” 

This  spreadsheet  is  very  similar  to  “P_Zeq_l.xls”  discussed  above.  It  allows  the  display  of 
the  surface  value  of  overpressure,  p(J,z  =  0) .  To  use  this  spreadsheet,  paste  the  three- 
column  output  in  “p2_zeq0.txt”  (Step  3)  into  columns  B,  C,  and  D.  Paste  the  two-column 
output  in  “p0_zeq0.txt”  (Step  3)  into  columns  F  and  G.  At  the  top  of  the  spreadsheet,  it  is 
also  necessary  to  enter  the  values  of  the  parameters:  M and  L. 

Calculations  proceed  exactly  as  for  the  above  spreadsheet,  with  the  final  result  of  total 
surface  overpressure,  p(T,z  =  0) ,  being  tabulated  in  column  L. 

4.3  Code  Inputs:  Summary  and  Utilization 

The  following  table  summarizes  the  inputs  needed  to  specify  a  particular  case. 


Input  No. 

Variable 

Units 

Input  in  Step(s): 

Input  Value  for 
Sample  Problem 

V.l 

Nx 

— 

1  and  5 

200 

V.2 

Xbegin 

— 

1  and  5 

-20.0 

V.3 

Xend 

— 

1  and  5 

20.0 

V.4 

M, 

— 

1  and  2 

1.5 

V.5 

z 

— 

1  and  5 

0.5 

V.6 

k 

— 

2 

14.0 

V.7 

A 

deg 

2 

0.0 

V.8 

T 

deg 

2 

0.0 

V.9 

Xstop 

— 

2  and  4 

20.0 

V.IO 

Xstep 

— 

2 

0.01 

V.ll 

Xi  begin 

— 

4 

-20.0 

V.12 

Xiend 

— 

4 

20.0 

V.13 

dxi 

— 

4 

0.01 

V.14 

P 

max 

psf 

Spreadsheets  2  and  3 

1.0 

V.15 

L 

ft 

Spreadsheets  2  and  3 

668. 

V.16 

S 

— 

Spreadsheets  2  and  3 

0.0267 

V.17 

ft/sec 

Assumed  =  1118.2 
Spreadsheets  2  and  3 

1118.2 

V.18 

Uoix) 

— 

1,  2,  3,  and  4 

N-wave 
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Inputs  V.4  -  V.8  and  V.14  -  V.17  are  defined  in  Section  2.1. 

Utilization  of  V.l  -  V.3  in  Steps  1  and  5 

The  first  three  inputs  (Nx,  Xbegin,  and  Xend)  allow  the  user  to  choose  the  x-interval  over 
which  solution  output  is  desired,  and  also  the  step  size  for  that  output.  For  example,  using 
the  Sample  Problem  values  of  (200,  -20.0,  and  20.0)  gives  output  from  x  =  -20  to  x  =  20  at 
intervals  of  A  x  =  40  /  200  =  0.2.  Note  that  these  three  inputs  are  not  used  for  any  numerical 
evaluations  of  integrals  or  derivatives. 

Evaluation  of  the  Integral  in  Eg.  (3.3)  in  Step  1 

The  integral  in  this  equation  is  evaluated  in  subroutine  s2d_genl,  which  is  contained  in  file 
“2d_gen26.for.”  [This  is  the  code  for  Step  1,  which  computes  the  Sawyers  (flat-surface) 
term.]  The  integral  is  evaluated  numerically  (see  next  section.)  The  code  employs  no  user 
inputs  to  establish  the  first-order  method  here.  The  interval  (0,1)  is  always  divided  evenly 
into  2000  intervals. 

N  =  2000  and  =  1/2000 
Numerical  Evaluations  in  Step  2 

The  inputs,  Xstop  (V.9)  and  Xstep  (V.IO),  are  used  in  the  routine  (Step  2)  which  computes 

A  A 

(p2{.x) .  Values  of  are  computed  starting  at  x  =  0  and  stopping  at  x  =  Xstop,  at  intervals  of 
Ax  =  Xstep.  Recommended  input  values  (see  Ref.  2)  are  Xstop  =  20  and  Xstep  =  0.01. 


A 

The  integral  S{x)  defined  in  Eqs.  (3.8)  and  (3.9)  is  evaluated  numerically  (see  next  section) 
with  an  integration  step  size  of  AXj  =  0.0001.  The  derivative  of  the  unit  surface  waveform 
which  appears  in  Eq.  (3.10)  is  computed  using  a  central  finite  difference  with  a  “grain”  of 
Ax,  =0.0001. 

Numerical  Evaluation  in  Step  3 

Values  of  p^.surf  [see  Eq.  (3.12)]  are  evaluated  in  Step  3  by  finite  differencing  consecutive 

values  of  which  were  computed  in  Step  2.  Thus  the  “grain”  of  the  finite  difference 
derivative  is  Ax  =  Xstep.  The  code  writes  only  every  fifth  value  to  output.  Thus,  the  usual 
situation,  after  execution  of  Step  3,  is  to  have  output  values  of  P2,surf  intervals  of  Ax  = 

0.05.  (NOTE:  The  calculation  of  the  surface  values  of  overpressure  due  to  wavy  interaction 
is  a  side  calculation.  The  results  are  not  used  in  subsequent  steps.) 

Numerical  Evaluations  in  Step  4 

The  variables  V.ll  -  V.13  are  used  in  the  routine  (Step  4)  which  computes  A(^) . 
Recommended  values  (see  Ref.  2)  are  Xibegin  =  -20,  Xiend  =  20,  and  dxi  =  0.01.  Values  of 
A  are  computed  starting  at  ^  =  Xibegin  and  stopping  at  ^  =  Xiend,  at  intervals  of 
A^  =  dxi . 
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The  integral  denoted  by  in  Eq.  (3.14)  is  computed  numerically  (see  next  section)  with 
A  X  =  Xstep.  For  numerical  evaluation,  the  upper  limit  on  the  integral 
(  X  =  oo )  is  taken  as  x  =  Xstop. 

The  integral  denoted  by  in  Eq.  (3.15)  is  computed  numerically  (see  next  section)  with 
A  X  =  0.01.  For  numerical  evaluation,  the  lower  and  upper  infinite  limits  on  the  integral  are 
taken  as  x  =  -Xstop  and  x  =  Xstop. 

The  function  q{x)  defined  in  Eq.  (3.16)  is  evaluated  using  a  forward  finite  difference  with 
Az  =  0.0001. 

Numerical  Evaluation  in  Step  5 

The  inverse  Fourier  transform  in  Eq.  (3.17)  is  evaluated  numerically  (see  next  section.)  For 
numerical  purposes,  the  integration  step  size  and  the  infinite  limits  on  the  integral  are  defined 
by  the  inputs  as  follows. 


A^=dxi,  =  Xibegin,  and  =  Xiend 

4.4  Numerical  Evaluation  of  Integrals 

There  are  five  integrals  evaluated  by  the  code.  These  are  given  in  Eqs.  (3.3),  (3.8),  (3.9), 
(3.14),  (3.15),  and  (3.17).  In  each  case,  the  following  method  is  used.  Suppose  the  integral 
of  interest  is: 


jfix)dx  (4.1) 

^u> 

If  the  integral  is  computed  using  N  evaluations  of  the  integrand,  then  define 

Ax  =  {Xf^,-Xijy)IN  (4.2) 

Xj=x^^+Ax/2+U-l)^  j  =  l,2,...,N  (4.3) 

fj  =  fixj)  j  =  (4.4) 

The  numerical  approximation  of  the  integral  is  then  given  by: 

N 

=  (4.5) 


For  the  five  integrals  of  interest,  the  values  of  the  integrand,  fj ,  are  complex  numbers. 
However,  the  variable  of  integration  (  denoted  byx,  x,  or  ^  )  is  always  real. 


20 


5  Code  Results  for  a  Sample  Problem 


A  sample  problem  is  now  presented  to  serve  as  a  test  case  for  the  code.  [To  reiterate  the 
statement  made  at  the  beginning  of  Section  4.2,  all  numerical  results  were  generated  using 
compiled  source  code.]  Suppose  a  sonic  boom  sweeps  over  an  area  of  deep  ocean  where 
the  wave  height  (trough-to-crest)  is  16  feet  and  the  wavelength  is  300  feet.  That  is, 

Ao=8  ft  (5.1) 

1  =  300  ft  (5.2) 

Suppose  also  that  the  boom  moves  along  the  water  surface  with  a  velocity  of  1677  ft/sec  in 
a  direction  exactly  opposite  to  that  of  the  propagation  of  the  water  waves.  That  is, 

U  =1611  ft/sec  (5.3) 


y  =  0“  (5.4) 

Assume  further  that  the  boom  signature  at  the  water  surface  is  an  N-wave,  668  feet  in 
length,  with  a  maximum  overpressure  of  1  psf,  and  that  we  are  interested  in  the  transmitted 


sound  directly  below  the  flight  path  of  the  supersonic  body. 

Uq  (x)  =  N-wave  (5.5) 

L  =  668  ft  (5.6) 

(5.7) 

A  =  0°  (5.8) 

From  the  above  parameters,  one  calculates  the  following  inputs  needed  for  running  the  code. 

=f//a^  =1677/1118.2  =  1.5  (5.9) 

k  =  {2^ /A)  1  =  2^-668/300  =  14  (5.10) 

S  =  AJ  A  =  8/300  =  0.0261  (5.11) 


Equations  (5.4)  through  (5.11)  specify  8  of  the  18  inputs  used  for  the  sample  problem. 
Values  used  for  the  remaining  10  inputs  are  shown  in  the  table  in  Section  4.3.  (For 
example:  Nx  =  200,  Xbegin  =  -20,  Xend  =  20,  Xstop  =  20,  Xstep  =  0.01,  dxi  =  0.01.) 
The  problem  was  run  for  four  different  scaled  depths:  z  =  0.01,  0.15,  0.5,  and  1.0 . 
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Output  results  for  the  flat  surface  (Sawyers)  term,  ,  are  shown  in  Fig.  5.  Results  for  the 
wavy  surface  term,  ,  are  shown  in  Fig.  6.  Results  for  the  total  disturbance,  Pi  +  P2,  are 
shown  in  Fig.  7. 


Time  (sec) 

Figure  6.  Wavy  surface  term,  pj  >  various  depths 


Overpressure  (psf)  Overpressure  (psf) 


Time  (sec) 

Figure  8.  Overpressures  at  z  =  1 . 


Time  (sec) 

Figure  9.  Overpressures  at  z  =  0.5 . 


Time  (sec) 

Figure  12.  Overpressures  at  z  =  0. 


In  addition  to  this  particular  sample  problem,  the  reader  may  consult  Ref.  5  for  additional 
sample  solutions. 


6  Summary 


The  sonic  boom  wavy  surface  code  (Ref.  1)  has  been  documented.  Explanations  have  been 
given  of  the  function  and  execution  of  the  various  code  modules.  Results  from  a  sample 
problem  have  been  plotted  in  order  to  provide  a  “benchmark”  for  future  users  and  against 
which  to  compare  possible  future  releases  of  this  code  by  Cheng  and  Lee. 
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